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The  SOR  iteratiqn  for  solving  linear  systems  of  equations  depends  upqn  an  overrelax¬ 
ation  factor  a.  We  show  that  for  the  standard  model  problem  of  Poisson js  equation  on  a 
rectangle,  the  optimal  u>  and  corresponding  convergence  rate  can  be  rigorously  obtained  by 
Fourier  analysis.  The  trick  is  to  tilt  the  space-time  grid  so  that  the  SOiy stencil  becomes 
symmetrical.  The  tilted  grid  also  gives  insight  into  the  relation  between  convergence  rates 
of  several  variants.  -N . .  .  , 
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1.  Introduction 


Fourier  analysis  has  been  used  for  nearly  fifty  years  to  teat  the  stability  of  time- 
dependent  finite  difference  formulas — the  “von  Neumann  method”  (9].  More  recently  it 
has  also  become  a  standard  tool  for  estimating  the  convergence  rate  of  multigrid  iterations 
[3].  But  the  analysis  of  the  more  classical  iteration  known  as  successive  overrelaxation — 
SOR — has  been  carried  out  by  other  means  [5,6,10,11].  The  reason  is  that  the  behavior 
of  SOR,  unlike  the  other  two  problems,  is  dominated  by  low-frequency  modes  that  are 
controlled  by  boundary  conditions.  The  obvious  application  of  Fourier  analysis  treats 
the  boundary  conditions  incorrectly,  and  leads  to  an  incorrect  prediction  of  the  optimal 
convergence  rate. 

In  this  note  we  show  that  if  the  computational  grid  is  tilted  by  a  certain  angle  in 
space  and  time,  Fourier  analysis  becomes  exact  for  the  standard  SOR  model  problem: 
the  five-point  discretization  of  Poisson’s  equation  on  a  rectangle  with  Dirichlet  boundary 
conditions,  with  the  variables  taken  in  the  natural  (typographical)  ordering. 

The  SOR  model  problem  was  first  analyzed  by  Frankel  in  1950  [6].  Our  approach 
leads  to  no  quantitative  results  that  Frankel  did  not  have,  but  makes  it  clear  why  the 
eigenvectors  of  the  SOR  iteration  have  the  form  they  do.  This  analysis  is  restricted  to 
the  rectangular  model  problem,  so  it  in  no  way  supplants  the  much  more  general  theory 
of  matrix  iterations  developed  by  ouag  in  the  early  1950’s  [10,11]. 

In  1956  Garabedian  proposed  a  t  analysis  of  SOR,  valid  in  the  limit  of  infinitesi¬ 
mal  mesh  spacing  [7].  He  observed  that  the  SOR  iteration  is  a  consistent  finite  difference 
approximation  of  a  time-dependent  partial  differential  equation,  so  that  its  rate  of  conver¬ 
gence  should  approximate  the  rate  of  decay  of  solutions  to  that  equation.  To  determine 
this  rate,  he  introduced  a  new  timelike  variable  a  =  t  +  z/2  +  y/2,  which  reduces  the 
differential  equation  to  a  canonical  form  that  can  be  analyzed  by  Fourier  methods.  Our 


tilting  of  the  grid  corresponds  exactly  to  Garabedian’s  introduction  of  the  variable  s. 
Thus  for  the  SOR  model  problem,  the  consideration  of  a  partial  differential  equation  is 
unnecessary,  and  indeed  the  analysis  in  the  discrete  domain  has  the  advantage  that  it  is 
exact  rather  than  approximate.  Garabedian’s  idea,  however,  provides  additional  insight 
and  is  applicable  to  more  general  problems. 

Approximate  Fourier  analysis  of  SOR  (on  the  usual  grid)  has  been  discussed  previ¬ 
ously  by  Kuo  [8]  and  Chan  and  Elman  [4]  and  probably  otheja^Our  tilted  grid  is  also 
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equivalent  to  the  “data  flow  times”  considered  by  Adams  and  Jordan  for  reasons  of  par* 
alienability  [1].  We  thank  all  of  these  authors,  and  also  David  Young,  for  their  advice 
and  encouragement. 


2.  Jacobi 

Consider  the  discrete  Poisson  problem 

+  uJ.k-l  +  “j.fc+l  _  4uy*)  =  fjkt  l<j,k<N-  1, 

Ujk  =  Fik,  j  =  0,  N  or  k  =  0,  N 

on  the  square  [O,*]2,  with  h  =  w/N.  Let  denote  the  approximation  to  the  exact 
discrete  solution  u  of  (1)  at  the  nth  step  of  an  iteration,  with  corresponding  error  = 

-  ujk- 

The  Jacobi  iteration  is  an  example  in  which  Fourier  analysis  works  straightforwardly 
[10,11],  The  errors  evolve  according  to 

“ft*  “*<*?->•*  +  +  SV.  +  <«.«)• 

t>;4  =  0,  i  =  0,AT  or  k  =  0,N. 

Let  us  consider  what  solutions  of  the  form  v"4  =  g(£,  rj)ne*^I+,,v)  this  iteration  admits  if 
we  ignore  the  boundary  conditions.  We  obtain  immediately 

9(Z,V)  =  $(«“**  +  *iih  +  «~ink  +  «inh)  =  |  (cos  £h  +  coer jh).  (3) 

This  is  the  amplification  factor  function  for  the  Jacobi  iteration.  The  essential  property 
is  that  it  is  an  even  function  of  (  and  rj: 

$(£.*))  =  =  $(-£,-»?)•  (4) 

If  we  take  as  initial  data  the  linear  combination 

sin£*sini7y  =  -i  («<(«*+•»»)  _  e«(-f*+w)  _  (5) 

where  £  and  r\  are  integers  in  the  range  l<£,r?<AT-l,  then  the  homogeneous  boundary 
conditions  are  satisfied  at  n  =  0.  By  (4),  it  follows  that  g(£,  rj)nBin£x  sinrjy  satisfies  both 
the  interior  formula  and  the  boundary  conditions  for  all  n  >  0,  and  therefore  sin£x  sinqy 
is  an  eigenvector  of  the  Jacobi  iteration  with  eigenvalue  ?((,*?).  Since  there  are  ( N  -  l)2 
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of  these  functions,  and  they  are  linearly  independent,  they  consitute  a  basis  for  the  set 
of  all  grid  functions  {v;*}.  Therefore  the  asymptotic  convergence  factor  for  the  Jacobi 
iteration  is  exactly 

J»cobi  _  max  |s(cos£h  +  cosn/i)!. 

l<f  ,rj  <N-l'2  n 

The  maximum  is  attained  with  £,r?  =  ±1  or  £,r)  =  ±[N  -  1), 


p Jacobi  _  cog/(l  ft,  1  -  Ifi2. 


(6) 


3.  SOR  and  Gauss-Seidel 

If  we  attempt  the  same  analysis  for  the  SOR  iteration, 

,n+l  —  (\  _  ,  .l.  —(..n+l  _i_  .,*»  _i_  „n+l 


«?:'  =  a  -  <-)«& + 


(7) 


the  result  is 


«(€,»?)  =  (!-«)  +  +  e-‘n*)  +  +  e**), 


that  is, 


*(£.»?) 


1  -  u  +  +  ei*k) 

1  -  |(e-‘^  +  «-••>*)  ‘ 


(8) 


Now,  (4)  no  longer  holds.  Therefore  $(£,r?)  does  not  give  us  the  eigenvalues  of  the  SOR 
iteration.  If  we  find  £  and  rj  to  maximize  \g(£ ,  »?)|,  there  is  no  reason  to  expect  the  resulting 
number  to  describe  the  convergence  of  SOR.  As  it  turns  out,  this  approach  produces  the 
correct  optimal  w,  to  leading  order  in  h,  but  a  convergence  rate  that  is  too  pessimistic  by 
a  factor  of  four  [4]. 

Figure  1  shows  how  the  situation  can  be  rescued.  Think  of  the  SOR  iteration  as 
inhabiting  a  regular  grid  in  two  space  and  one  time  dimensions  ( j,k,n ).  Its  stencil  con¬ 
nects  six  points  in  an  asymmetrical  pattern,  or  four  points  in  the  one-dimensional  case 
portrayed  in  the  figure.  Because  of  this  asymmetry,  (4)  does  not  hold.  But  if  we  introduce 
the  new  “time”  index 

u  =  2n  +  j  +  k,  (9) 


the  stencil  becomes  symmetrical.  Let  us  look  for  solutions  to  (7)  of  the  form 


v,k 


=  y(e,r?)V({*+,"')- 


(10) 
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Figure  1.  The  SOR  stencil  superimposed  on  a  space-time  grid  (one  space  di¬ 
mension).  Introduction  of  the  “tilted”  time  index  v  makes  the  stencil  symmetric, 
so  that  Fourier  analysis  can  be  applied.  The  red-black  labels  are  explained  in 
Section  4. 


In  the  j,k,u  variables  (7)  becomes 

<#’  =  (i  -  "K* + }«.■.. + tm 

and  so  a  suitable  value  for  y(£,  rj)  is  either  root  of  the  quadratic  equation 

=  (1  -  <*/)  +  ^(cos$h  +  cosqh)ff($,»7).  (12) 

For  each  £  and  rj  we  now  have  a  pair  of  amplification  factors  g±(^,rj),  and  they  satisfy 
the  symmetry  condition  (4).  Therefore  for  any  integers  rj  in  the  range  1  <  rj  <  N  - 1, 
the  functions  g(Z,n)y sinfxsin»7y  are  eigenmodes  of  the  SOR  iteration  in  the  u  direction. 
To  speak  in  terms  of  eigenvectors,  we  note  that  SOR  is  a  two-step  formula  with  respect  to 
v,  but  it  can  be  recast  as  a  one-step  iteration  («*'“*,  v1'-1)  (vl',v*'+1)  with  eigenvectors 

(sin^zsinqy ,  $(£,rj)sin£xsinr7y)  and  eigenvalues  y(^,f?)s. 

It  takes  two  steps  in  v  to  advance  one  step  in  n.  We  conclude  that  the  asymptotic 
convergence  factor  for  SOR  is  exactly 

^OR  =  ,*/"25r  .  W  !**(*»  *)!*■  (13) 

In  the  original  j,k,n  coordinates,  the  eigenmodes  become 

v"*  =  9  (£,  q)ln+>+‘ginfz  sinijy, 
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and  the  corresponding  eigenvectors  are 

g{S,riy+isinZxainT)y. 

This  matches  the  results  of  Frankel  and  others  derived  by  different  means. 

All  that  remains  is  algebra.  For  any  £,  rj,  and  ui  >  1,  the  solutions  to  (12)  are 

j  fjJ 

9±\Z,V)  =  a±  yj a2  -  (w  -  1),  a=  -(cosfh  +  cosrjh), 
and  the  larger  of  these  two  numbers  has  magnitude 


max\g±{Z,ri)\  = 
'  > 


|a|  +  y/a1  -  (<jj  -  1)  if  a2  >  w  -  1, 


/ui  —  1 


if  a2  <  w  —  1  . 


For  fixed  w,  this  quantity  can  evidently  be  maximized  with  respect  to  (  and  17  by  taking 
£  =  17  =  1  (among  other  values)  and  a  =  ycos/i. 

The  Gauss- Seidel  iteration  corresponds  to  u  =  1.  In  this  case  (16)  becomes  2)a)  = 
cos/i,  so  by  (13)  we  have 

pGS  =  cos2h  w  1  -  A2.  (17) 

To  find  the  optimal  overrelaxation  factor  for  SOR,  we  examine  the  dependence  of  (16) 
on  w  with  £  =  fj  =  1.  It  is  readily  verified  that  values  u  >  2  lead  to  p®°R  >  ®o  they  are 
out  of  the  running,  and  we  can  assume  1  <  w  <  2.  In  this  range  the  second  line  of  (16) 
obviously  increases  with  w,  and  differentiation  confirms  that  the  first  line  decreases  with 
<jj.  Therefore  the  optimal  w  is  the  crossover  value  w  -  1  =  a2  =  (^co ah)7,  which  reduces 
after  a  little  work  to 

=  TTm  ca2'2h-  (I8) 

By  (13)  and  (16),  the  corresponding  convergence  rate  is 

Pop.  -  wopt  -  1  -  l+ainhPal  2h- 


(19) 


4.  Relating  various  methods 


The  change  to  tilted  coordinates  has  the  additional  advantage  of  clarifying  the  re¬ 
lationships  between  convergence  rates  of  different  iterative  methods.  For  example,  it  is 
well  known  that  the  Gauss-Seidel  iteration  is  twice  as  fast  as  Jacobi,  as  is  confirmed  by 
comparing  (6)  and  (17).  The  tilted  coordinates  provide  a  simple  explanation  of  why  this 
is  so.  Gauss-Seidel  corresponds  to  the  case  w  =  1  of  (11),  and  this  is  precisely  the  Jacobi 
iteration  with  respect  to  v.  The  factor  of  two  comes  from  the  fact  that  it  takes  two  steps 
in  v  to  advance  one  step  in  n. 

As  another  example,  consider  the  SOR  iteration  with  the  grid  points  taken  in  the 
red-black  or  checkerboard  ordering.  This  means  that  the  v,*  with  j  +  k  even  (red  points) 
are  processed  before  the  Vjt  with  j +  k  odd  (black  points),  and  the  iteration  takes  the 
form 


-  (!  ~  u)vjk  +  j(v?~l,k  +  u?+l,*  +  v?,k-i  +  v?,k+i) 
=  (1  -  «)«&  +  + «£&) 


for  j  +  k  even, 


for  j  +  k  odd. 


For  each  u>,  this  method  has  the  same  convergence  rate  as  the  iteration  (7)  in  the  natural 
ordering.  Young  proved  this  algebraically  by  determining  the  eigenvalues  and  eigenvectors 
of  the  associated  iteration  matrices  [11].  Again,  the  change  to  tilted  coordinates  gives  a 
more  intuitive  explanation,  as  illustrated  in  Figure  1  for  the  case  of  one  space  dimension. 
At  step  v  we  are  computing  only  the  red  points,  say,  and  at  step  t/+l  only  the  black  points. 
Recasting  SOR  as  a  one-step  iteration  (v1'-1, t/*'-1)  *-*  (v", ul/+1),  as  in  the  last  section, 
we  obtain  simply  the  red-black  ordering.  Thus  Figure  1  can  be  viewed  as  depicting  an 
SOR  iteration  either  in  j,n  coordinates  with  the  natural  ordering,  or  in  j,v  coordinates 
with  the  red-black  ordering.  Hence  these  two  orderings  must  have  the  same  asymptotic 
convergence  rate. 


The  conclusions  above  depend  on  the  fact  that  the  convergence  rate  is  independent 
of  the  particular  initial  data  used,  depending  only  on  the  eigenvectors.  Note  that  we  can 
switch  back  and  forth  between  arbitrary  data  at  fixed  n  or  at  fixed  i>  by  taking  partial 
iterations  over  a  triangular  portion  of  the  grid.  In  fact,  writing  out  these  partial  iterations 
algebraically  gives  a  similarity  transformation  relating  the  iteration  matrices. 

In  this  paper  we  have  considered  just  the  five-point  Poisson  mode!  problem,  and 
presented  an  easy  way  to  obtain  classical  results  with,  we  hope,  additional  insight.  The 
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tilted  grid  may  also  prove  useful  in  obtaining  new  results.  It  has  already  been  applied  to 
settle  a  conjecture  of  Adams  and  Jordan  regarding  the  equivalence  of  certain  orderings 
for  the  nine-point  stencil  (l|.  These  results  will  be  reported  in  [2]. 


References 

[1]  L.  M.  Adams  and  H.  F.  Jordan,  Is  SOR  color-blind?,  SIAM  J.  Sci.  Sfcat.  Comp.  7 
(1986),  490-506. 

[2]  L.  M.  Adams,  R.  J.  LeVeque,  and  D.  M.  Young,  Analysis  of  the  SOR  iteration  for  the 
9-point  Laplacian,  to  appear  as  an  ICASE  report. 

[3]  A.  Brandt,  Multi-level  adaptive  solutions  to  boundary-value  problems,  Math.  Comp. 
31  (1977),  333-390. 

[4]  T.  F.  Chan  and  H.  C.  Elman,  Fourier  analysis  of  preconditioned  iterative  methods, 
Yale  Univ.  Dept,  of  Comp.  Sci.  Rep.  Yaleu/DCS/RR-410,  to  appear. 

[5]  G.  E.  Forsythe  and  W.  R.  Wasow,  Finite- Difference  Methods  for  Partin]  Differential 
Equations,  Wiley,  1960. 

[6]  S.  Frankel,  Convergence  rates  of  iterative  treatments  of  partial  differential  equations, 
Math.  Comp.  4  (1950),  65-75. 

[7]  P.  Garabedian,  Estimation  of  the  relaxation  factor  for  small  mesh  size,  Math.  Comp. 

10  (1956),  183-185. 

[8]  C.-C.  Kuo,  B.  C.  Levy,  and  B  R.  Musicus,  A  local  relaxation  method  for  solving 
elliptic  PDEs  on  mesh-connected  processor  arrays,  SIAM  J.  Sci.  Stat.  Comp.,  to 
appear. 

[9]  R.  D.  Richtmyer  and  K.  W.  Morton,  Difference  Methods  for  Initial-Value  Problems, 
2nd  ed.,  Wiley-Interscience,  1967. 

[10]  R.  S.  Varga,  Matrix  Iterative  Analysis,  Prentice-Hall,  1962. 

[11]  D.  Young,  Iterative  Solution  of  Large  Linear  Systems,  Academic  Press,  1971. 


7 


Standard  Bibliographic  Page 


1  Report  No  NASA  CR- 178191  i  2.  Government  Acceeeion  No. 

ICASE  Report  No.  86-63 _ j _ 

4  Title  and  Subtitle 

FOURIER  ANALYSIS  OF  THE  SOR  ITERATION 


7  Author(s) 

Randall  J.  LeVeque,  Lloyd  N.  Trefethen 


9  Performing  Organization  Name  and  Address 

Institute  for  Computer  Applications  in  Science 
and  Engineering 

Mail  Stop  132C,  NASA  Langley  Research  Center 
Hampton r  VA  23665-522,5 _ 

12.  Sponsoring  Agency  Name  and  Address 


National  Aeronautics  and  Space  Administration 
Washington.  D.C.  20546 _ 

15.  Supplementary  Notes 


,  3.  Recipient  *  Catalog  No. 

!  ! 

■  —  -  '  —  - -j 

i  5.  Report  Date 

t 

! _ Sp.ptemher  1986 _ I 

I  6.  Performing  Organization  Code 

'  i 

4 _ J 

8.  Performing  Organization  Report  No.  ! 

I _ 86-63 _ I 

4  10.  Work  Unit  No 

i 

11.  Contract  or  Grant  No. 

NAS1-181Q7 _ ' 

J  13.  Type  of  Report  and  Penod  Covered 

— Coat r actor  Report - 

14.  Sponsoring  Agency  Code 

- — 303-90-21—01-  ■  -  -  - 


Langley  Technical  Monitor:  SIAM  Num.  Analysis 

J.  C.  South 


Final  Report 

16  Abstract 


The  SOR  iteration  for  solving  linear  systems  of  equations  depends  upon  an  j 
overrelaxation  factor  1  .  We  show  that  for  the  standard  model  problem  of  j 
Poisson's  equation  on  a  rectangle,  the  optimal  1  and  corresponding  convergence 
rate  can  be  rigorously  obtained  by  Fourier  analysis.  The  trick  is  to  tilt  the  j 
space-time  grid  so  that  the  SOR  stencil  becomes  symmetrical.  The  tilted  grid 
also  gives  Insight  Into  the  relation  between  convergence  rates  of  several 
variants . 
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